import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

/**
 * 使用前向算法计算T时刻观测到某一观测序列的概率
 * 前向算法是求T时刻每种隐藏状态观测到指定观测序列的概率之和
 * 利用动态规划，从1时刻往后推，1时刻隐藏状态为i观测序列为y1的概率为初始隐藏状态概率πi乘以生成观测状态概率Biy1，
 * 那么t+1时刻隐藏状态为i观测序列为y1…yn的概率为t时刻所有状态的观测到y1…yn-1然后向i状态转移并观测到yn的概率之和，
 * 最终p(O|λ)就等于T时刻所有可能的隐藏状态观测到观测序列为y1…yn的概率之和
 */
public class ForwardAlgorithm {

    static int[] Q = {1, 2, 3}; //隐藏状态集合，

    static List<Integer> V = Arrays.asList(1, 2); //观测状态集合

    static double[] PAI = {0.2, 0.4, 0.4}; //初始隐藏状态概率

    static double[][] A = {{0.5, 0.2, 0.3}, {0.3, 0.5, 0.2}, {0.2, 0.3, 0.5}}; //隐藏状态的状态转移矩阵

    static double[][] B = {{0.5, 0.5}, {0.4, 0.6}, {0.7, 0.3}}; //生成观测状态概率矩阵

    static int T = 3; //T时刻

    static int[] O = {1, 2, 1}; //T时刻得到的观测序列

    public static void main(String[] args) {
        System.out.println("时刻"+T+"观测序列为"+Arrays.toString(O)+"的概率为：" + calculateProbability());
        long experimentTimes = 100000;
        System.out.println("实验"+experimentTimes+"次，观测序列为"+Arrays.toString(O)+"的比率为：" + extractExperiment(100000));
    }

    /**
     * 使用前向算法计算T时刻观测到观测序列O的概率
     */
    private static double calculateProbability() {
        //计算T=1时刻每种隐藏状态观测到O[0]的概率
        int o1Index = V.indexOf(O[0]);
        List<Double> pList = new ArrayList<>();
        for (int i = 0;i < Q.length;i++) {
            pList.add(PAI[i] * B[i][o1Index]);
        }

        for (int t=2;t <= T;t++) {
            List<Double> currentPList = new ArrayList<>();
            for (int i = 0;i < Q.length;i++) {
                //计算t时刻向隐藏状态Q[i]转移的概率tp
                double tp = 0;
                for (int j = 0;j < Q.length;j++) {
                    tp += pList.get(j)*A[j][i];
                }
                int oIndex = V.indexOf(O[t-1]);
                //t时刻隐藏状态为Q[i]并且观测序列为O[0]到O[t-1]的概率
                currentPList.add(tp*B[i][oIndex]);
            }
            pList = currentPList;
        }

        double result = 0;
        for (int i = 0;i < pList.size();i++) {
            result += pList.get(i);
        }
        return result;
    }

    /**
     * 模拟真实情况得到的给定观测序列O的比率
     * @param experimentTimes 实验次数
     */
    private static double extractExperiment(int experimentTimes) {
        int hitCounts = 0;
        for (int i=0;i<experimentTimes;i++) {
            List<Integer> resultList = new ArrayList<>();
            int qIndex = randomWithPs(PAI);
            int q1 = Q[qIndex]; //初始隐藏状态
            int v1 = V.get(randomWithPs(B[qIndex])); //初始生成观测状态
            resultList.add(v1);
            for (int t = 2;t<=T;t++) {
                qIndex = randomWithPs(A[qIndex]);
                resultList.add(V.get(randomWithPs(B[qIndex])));
            }
            boolean equal = true;
            for (int j = 0;j < resultList.size();j++) {
                if (O[j] != resultList.get(j)) {
                    equal = false;
                }
            }
            if (equal) {
                hitCounts++;
            }
        }
        return 1.0*hitCounts/experimentTimes;
    }

    /**
     * 在一个概率数组中，以下标所在概率数值为随机权重，随机返回一个下标
     * @param ps
     * @return
     */
    private static int randomWithPs(double[] ps) {
        double rNumber = Math.random();
        double p = 0;
        for (int i=0;i<ps.length;i++) {
            p = p + ps[i];
            if (rNumber < p) {
                return i;
            }
        }
        return -1;
    }
}
